{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Imports\n",
    "\n",
    "# utility modules\n",
    "import glob\n",
    "import os\n",
    "import sys\n",
    "import re\n",
    "\n",
    "# the usual suspects:\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# specialty modules\n",
    "import h5py\n",
    "import pyproj\n",
    "\n",
    "# modules you'll need if you're downloading the data:\n",
    "from icepyx import icesat2data as ipd\n",
    "import shutil\n",
    "import geopandas as gpd\n",
    "\n",
    "# run matplotlib in 'widget' mode\n",
    "%matplotlib widget\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# # Jessica's function for loading multiple data requests\n",
    "# #define function to loop through series of requests to get data\n",
    "\n",
    "# def multiple_is2_requests(earthdata_email, earthdata_uid, data_home, requests, subset, download):\n",
    "#     for req in requests:\n",
    "#         #### Download this data: (uncomment and run)\n",
    "#         region_a = ipd.Icesat2Data(req['short_name'], req['spatial_extent'], req['date_range'])\n",
    "#         if download==True:\n",
    "#             region_a.earthdata_login(earthdata_uid, earthdata_email)\n",
    "#             region_a.download_granules(data_home,subset=subset)\n",
    "#         print(region_a.dataset)\n",
    "#         print(region_a.dates)\n",
    "#         print(region_a.start_time)\n",
    "#         print(region_a.end_time)\n",
    "#         print(region_a.dataset_version)\n",
    "#         print(region_a.spatial_extent)\n",
    "#         region_a.visualize_spatial_extent()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A note on data\n",
    "Since the hackweek is over, you'll need to download some ATL03 and ATL06 yourself put it somewhere sensible.  Mine is in /home/jovyan/tutorial-data, but you may need to edit this cell to match where your data ended up.  The next cell contains the code to download the data from NSIDC.  It takes a while to run, so if you run the cell this may be a good time to get a cup of coffee.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# earthdata_email='your email here'\n",
    "# earthdata_name='your earthdata login'\n",
    "\n",
    "# data_home='/home/jovyan/tutorial-data'\n",
    "# #set DOWNLOAD to true if you're ready to download the data:\n",
    "# download=True\n",
    "# subset=True\n",
    "    \n",
    "# requests=[\n",
    "#     {   'short_name' : 'ATL06',\n",
    "#         'spatial_extent' :[-102, -76, -98, -74.5],\n",
    "#         'date_range' : ['2018-10-14','2020-04-01']},\n",
    "#     {   'short_name' :'ATL03',\n",
    "#         'spatial_extent' :[-102, -76, -98, -74.5],\n",
    "#         'date_range' : ['2018-10-14','2020-04-01']},\n",
    "#     {   'short_name' :'ATL03',\n",
    "#         'spatial_extent' :[-102, -76, -98, -74.5],\n",
    "#         'date_range' : ['2019-12-28', '2019-12-29']}]\n",
    "\n",
    "# multiple_is2_requests(earthdata_email, earthdata_uid, data_home, requests, subset, download);\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Land ice applications: Getting familiar with ICESat-2 products over land ice\n",
    "\n",
    "## 1. Introduction\n",
    "In this tutorial, we're going to take a first look at some of the ICESat-2 data over land ice.  We're going to demonstrate how ALT06 segments correspond to ATL03 photons, and demonstrate the repeat-track structure of ATL06 data. '\n",
    "\n",
    "## 1.1 Learning goals\n",
    "Our learning goals include understanding:\n",
    "* how ice-sheet surfaces are represented by ATL03 photons\n",
    "* how ATL06 elevations correpond to ATL03 photon clouds\n",
    "* how small-scale features appear in ATL06\n",
    "* the repeat structure of ATL06\n",
    "* how cross-track slope can affect elevation differences\n",
    "and we'll also give:\n",
    "* a sneak peak at the ATL11 product.\n",
    "\n",
    "## 1.2 Tools presented and developed:\n",
    "Along the way, we'll present:\n",
    "* A reader for ATL03\n",
    "* A reader for ATL06\n",
    "* Geographic projections with pyproj\n",
    "* Filtering ATL06 elevations using slope and quality parameters\n",
    "* Mapping surface elevations and profiles with matplotlib\n",
    "* Visualizing ATL03 and ATL06 together."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2.1 Modules Used:\n",
    "If you look at the first cell in the notebook, you can see that we're importing:\n",
    "* h5py: a library to read (and write) hdf5 data\n",
    "* pyproj: a library of geographic corrections\n",
    "\n",
    "We'll also us my pointCollection package, which contains code to read a variety of datatypes, which we'll use to read image data and to take a first look at the ATL11 product."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# try:\n",
    "#     import pointCollection as pc\n",
    "# except Exception:\n",
    "#     !python3 -m pip install --user git+https://github.com/smithb/pointCollection.git\n",
    "import pointCollection as pc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll also use Tyler Sutterley's excellent read_HDF5_ATL03 package (stored in the 'readers' directory, one level above the current directory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "sys.path.append(os.path.join(os.getcwd(), '..'))\n",
    "from readers.read_HDF5_ATL03 import read_HDF5_ATL03\n",
    "from readers.get_ATL03_x_atc import get_ATL03_x_atc\n",
    "data_root = os.getcwd()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2.2 Data used\n",
    "In this tutorial, we'll use the ATL03 and ATL06 data that were downloaded in the top couple of cells.  You'll have put them in the data_root directory specified in those cells, and we'll look for our mosaic background image in the same location."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Image background\n",
    "\n",
    "We'll also use an image mosaic as background for Antarcitca.  It's available in the shared space on jovyan in the MOA subdirectory\n",
    "\n",
    "If you're running the script outside the hackweek.io world, you'll need to download the files from NSIDC.  I've posted a 1-km downscale of the raw (125-meter) mosaic, but the 750-meter mosaic here will do nicely:\n",
    "\n",
    "https://daacdata.apps.nsidc.org/pub/DATASETS/nsidc0593_moa2009/geotiff/moa750_2009_hp1_v01.1.tif.gz\n",
    "\n",
    "Try downloading the linked file (you'll need to log in to Earthdata to get it), then put it in the data_root directory on your home machine."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.3 Geographic setting\n",
    "\n",
    "Let's read the the Mosaic of Antarctica for the region around Pine Island Glacier.  The mosaic is in polar-stereographic coordinates, so we'll need to project the geographic coordinates of the box into that projection to decide what part of the mosaic to read."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# spatial_extent = np.array([-102, -76, -98, -74.5])\n",
    "# lat=spatial_extent[[1, 3, 3, 1, 1]]\n",
    "# lon=spatial_extent[[2, 2, 0, 0, 2]]\n",
    "# print(lat)\n",
    "# print(lon)\n",
    "# # project the coordinates to Antarctic polar stereographic\n",
    "# xy=np.array(pyproj.Proj(3031)(lon, lat))\n",
    "# # get the bounds of the projected coordinates \n",
    "# XR=[np.nanmin(xy[0,:]), np.nanmax(xy[0,:])]\n",
    "# YR=[np.nanmin(xy[1,:]), np.nanmax(xy[1,:])]\n",
    "\n",
    "# ### EDIT THIS PATH TO THE MOA LOCATION IF YOU NEED TO\n",
    "# MOA=pc.grid.data().from_geotif(os.path.join(data_root, 'MOA','January_original.tif'), bounds=[XR, YR])\n",
    "\n",
    "# # show the mosaic:\n",
    "# plt.figure()\n",
    "# MOA.show(cmap='gray', clim=[14000, 17000])\n",
    "# plt.plot(xy[0,:], xy[1,:])\n",
    "# plt.title('Mosaic of Antarctica for Pine Island Glacier')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the Pine Island Glacier.  It flows from the top of the map to the bottom (actually east to west, but it's rotated in the Polar Stereographic projection used here).  The striped area running from the middle of the page to the bottom is the Pine Island Ice Shelf, where the ice is floating.  We expect to see crevasses at the edges of the ice shelf (or really, anywhere in the ice shelf).  This area is often cloudy, so we expect to have trouble seeing the surface a lot of the time.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. Data in along-track format\n",
    "\n",
    "# 2.1 ATL03 elevations\n",
    "\n",
    "Before we start looking at the ATL06 data we've downloaded, let's have a look at some of the ATL03 data that were used to make them.  One of the source ATL03 files is in the shared folder, and we'll read it with Tyler Sutterley's excellent \"read_HDF5_ATL03\" function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read the data:\n",
    "rgt=\"0027\"\n",
    "cycle=\"06\"\n",
    "# read the IS2 data with Tyler's ATL03 reader:\n",
    "sys.path.append(os.path.join(os.getcwd(), '..'))\n",
    "from readers.read_HDF5_ATL03 import read_HDF5_ATL03\n",
    "from readers.get_ATL03_x_atc import get_ATL03_x_atc\n",
    "\n",
    "\n",
    "ATL03_file='ATL03_20190128163325_04770206_003_01.h5'\n",
    "IS2_atl03_mds, IS2_atl03_attrs, IS2_atl03_beams =read_HDF5_ATL03(ATL03_file)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The data are returned in a set of dictionaries that mimic the structure of an ATL03 file.  To help visualize the data, we're going to calculate an along-track coordinate for every photon in the cloud (x_atc).  This is a slightly complex job, and there's a helper function in the readers directory that you can look at if you want the details."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# add x_atc to the ATL03 data structure (this function adds to the LS2_ATL03_mds dictionary)\n",
    "get_ATL03_x_atc(IS2_atl03_mds, IS2_atl03_attrs, IS2_atl03_beams)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1.1 Plotting ATL03 photons\n",
    "Now let's plot the ATL03 photons.  We'll plot all the photon heights as small black dots, then plot the photons that the ATL03 land-ice signal finder designates as surface (with low, medium, or high confidence) in green.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d96519a7c994463f8fecbcdd2d6590c7",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#-- select the 2l beam from ATL03\n",
    "D3 = IS2_atl03_mds['gt2l']\n",
    "\n",
    "#-- create scatter plot of photon data (e.g., photon elevation vs x_atc)\n",
    "fig=plt.figure(figsize=(6,4))\n",
    "ax = fig.add_subplot(111)\n",
    "ax.plot(D3['heights']['x_atc'], D3['heights']['h_ph'],'k.',markersize=0.25, label='all photons')\n",
    "LMH=D3['heights']['signal_conf_ph'][:,3] >= 2\n",
    "ax.plot(D3['heights']['x_atc'][LMH], D3['heights']['h_ph'][LMH],'g.',markersize=0.5, label='flagged photons')\n",
    "h_leg=ax.legend()\n",
    "\n",
    "ax.set_xlabel('x_atc, m')\n",
    "ax.set_ylabel('h, m')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What we see here is a fat black bar, representing all the photons that were close to the surface (in a ~200 m window), with a green line representing the photons flagged by ATL03 as representing the surface. If we zoom in (using the square button on the right) we can see that the black bar is made up of lots of individual photons, with a darker line that represents the surface.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.2 ATL06 elevations\n",
    "\n",
    "ATL03 can show where the surface is, but it doesn't explicitly report a surface height.  It's also a bulky product that can be slow to load.  ATL06 reduces elevation data to a 20-meter posting, and gives one elevation per 20 meters (instead of hundreds in ATL03). \n",
    "\n",
    "### 2.2.1 simple ATL06 reader\n",
    "To help work with ATL06 data, we'll use a simple piece of code that reads ATL06 data one beam at a time and stores it in a dictionary.  The code has a default set of fields to be read from an ATL06 file, and stores the data from each field in the output dictionary.  It also takes care of removing bad data, by setting values that are marked as invalid in the file to NaN.  It also uses hte \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def atl06_to_dict(filename, beam, field_dict=None, index=None, epsg=None):\n",
    "    \"\"\"\n",
    "        D6=atl06_to_dict(ATL06_file,'/gt2l', index=None, epsg=3031)\n",
    "\n",
    "        Read selected datasets from an ATL06 file\n",
    "\n",
    "        Input arguments:\n",
    "            filename: ATl06 file to read\n",
    "            beam: a string specifying which beam is to be read (ex: gt1l, gt1r, gt2l, etc)\n",
    "            field_dict: A dictinary describing the fields to be read\n",
    "                    keys give the group names to be read, \n",
    "                    entries are lists of datasets within the groups\n",
    "            index: which entries in each field to read\n",
    "            epsg: an EPSG code specifying a projection (see www.epsg.org).  Good choices are:\n",
    "                for Greenland, 3413 (polar stereographic projection, with Greenland along the Y axis)\n",
    "                for Antarctica, 3031 (polar stereographic projection, centered on the Pouth Pole)\n",
    "        Output argument:\n",
    "            D6: dictionary containing ATL06 data.  Each dataset in \n",
    "                dataset_dict has its own entry in D6.  Each dataset \n",
    "                in D6 contains a numpy array containing the \n",
    "                data\n",
    "    \"\"\"\n",
    "    if field_dict is None:\n",
    "        field_dict={None:['latitude','longitude','h_li', 'atl06_quality_summary'],\\\n",
    "                    'ground_track':['x_atc','y_atc'],\\\n",
    "                    'fit_statistics':['dh_fit_dx', 'dh_fit_dy']}\n",
    "    D={}\n",
    "    file_re=re.compile('ATL06_(?P<date>\\d+)_(?P<rgt>\\d\\d\\d\\d)(?P<cycle>\\d\\d)(?P<region>\\d\\d)_(?P<release>\\d\\d\\d)_(?P<version>\\d\\d).h5')\n",
    "    with h5py.File(filename,'r') as h5f:\n",
    "        for key in field_dict:\n",
    "            for ds in field_dict[key]:\n",
    "                if key is not None:\n",
    "                    ds_name=beam+'/land_ice_segments/'+key+'/'+ds\n",
    "                else:\n",
    "                    ds_name=beam+'/land_ice_segments/'+ds\n",
    "                if index is not None:\n",
    "                    D[ds]=np.array(h5f[ds_name][index])\n",
    "                else:\n",
    "                    D[ds]=np.array(h5f[ds_name])\n",
    "                if '_FillValue' in h5f[ds_name].attrs:\n",
    "                    bad_vals=D[ds]==h5f[ds_name].attrs['_FillValue']\n",
    "                    D[ds]=D[ds].astype(float)\n",
    "                    D[ds][bad_vals]=np.NaN\n",
    "    if epsg is not None:\n",
    "        xy=np.array(pyproj.proj.Proj(epsg)(D['longitude'], D['latitude']))\n",
    "        D['x']=xy[0,:].reshape(D['latitude'].shape)\n",
    "        D['y']=xy[1,:].reshape(D['latitude'].shape)\n",
    "    temp=file_re.search(filename)\n",
    "    D['rgt']=int(temp['rgt'])\n",
    "    D['cycle']=int(temp['cycle'])\n",
    "    D['beam']=beam\n",
    "    return D"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll use this to read the ATL06 data that correspond to the ATL03 data in the previous cell, and will plot the ATL06 elevations on top of the ATL03 photons."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# find the matching ATL06 file\n",
    "ATL06_file='ATL06_20190128163325_04770206_003_01.h5'\n",
    "D6=atl06_to_dict(ATL06_file,'/gt2l', index=None, epsg=3031)\n",
    "\n",
    "# plot the elevations on top of the previous axes.  You should be able to scroll up to the previous plot and see the ATL06 points.\n",
    "ax.plot(D6['x_atc'], D6['h_li'],'r.', label='ATL06')\n",
    "ax.legend();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'x_atc': 83928, 'h_li': 83928}"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "{field:D6[field].size for field in ['x_atc','h_li']}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1.2 Cloudy tracks and singnal-finding blunders"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cycle 6 happened in pretty nice weather.  Not so much cycle 5.  Let's take a look at the ATL03 and the ATL06 for cycle 5 and see what happens when things don't go well for ICESat-2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ATL03_20190128163325_04770206_003_01.h5\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "ab92a14d9ddf496f917a83f83be7a0e0",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "rgt=\"0027\"\n",
    "cycle=\"05\"\n",
    "beam='gt2r'\n",
    "\n",
    "# read the IS2 data with Tyler's ATL03 reader:\n",
    "ATL03_file='ATL03_20190128163325_04770206_003_01.h5'\n",
    "print(ATL03_file)\n",
    "IS2_atl03_mds, IS2_atl03_attrs, IS2_atl03_beams =read_HDF5_ATL03(ATL03_file)\n",
    "# add x_atc to the ATL03 data structure (this function adds to the LS2_ATL03_mds dictionary)\n",
    "get_ATL03_x_atc(IS2_atl03_mds, IS2_atl03_attrs, IS2_atl03_beams)\n",
    "\n",
    "#-- select the 2r beam from ATL03\n",
    "D3 = IS2_atl03_mds[beam]\n",
    "\n",
    "# find the matching ATL06 file\n",
    "ATL06_file='ATL06_20190128163325_04770206_003_01.h5'\n",
    "D6=atl06_to_dict(ATL06_file, beam, index=None, epsg=3031)\n",
    "#-- create scatter plot of photon data (e.g., photon elevation vs x_atc)\n",
    "%matplotlib widget\n",
    "f1,ax = plt.subplots(num=1,figsize=(6,4))\n",
    "TEP=(D3['heights']['signal_conf_ph'][:,3] <-1)   \n",
    "ax.plot(D3['heights']['x_atc'][TEP], D3['heights']['h_ph'][TEP],'b.',markersize=0.25, label='TEP')\n",
    "BG=(D3['heights']['signal_conf_ph'][:,3] ==0)   |  (D3['heights']['signal_conf_ph'][:,3] ==1)\n",
    "ax.plot(D3['heights']['x_atc'][BG], D3['heights']['h_ph'][BG],'k.',markersize=0.25, label='Background')\n",
    "LMH=D3['heights']['signal_conf_ph'][:,3] >= 2\n",
    "ax.plot(D3['heights']['x_atc'][LMH], D3['heights']['h_ph'][LMH],'g.',markersize=0.5, label='flagged photons')\n",
    "ax.plot(D6['x_atc'], D6['h_li'],'r.', label='ATL06')\n",
    "h_leg=ax.legend()\n",
    "\n",
    "plt.title(os.path.basename(ATL03_file))\n",
    "\n",
    "ax.set_xlabel('x_atc, m')\n",
    "ax.set_ylabel('h, m')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now the window of photons around the surface is much more scattered, there are some places where ATL06 is reporting heights for something that's probably not the surface, and there's a stripe of extra photons below the surface that doesn't make any obvious sense.  The new (blue) stripe of photons is the TEP or Transmitter Echo Pulse, which records a packet of photons that are looped from the transmit to the receive side of ATLAS to help monitor its impulse response.  ATL06 ignores these photons, and we should too.  \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1.3 Filtering out bad returns with ATL06_quality summary and with along-track differencing.\n",
    "\n",
    "ATL06 \"points\" aren't just points, though.  They contain a variety of different flags that can help us filter out bad data.  Two of those that we can use are:\n",
    "\n",
    "-- atl06_quality_summary : a flag that evaluates whether any problems have been found for a particular segment.  A 'zero' value indicates that the automatic filters haven't found anything wrong with the segment.\n",
    "-- dh_fit_dx : The surface slope estimated for the segment in the along-track direction.  We can use this to check whether each segment is consistent with its neighbors.\n",
    "\n",
    "Let's add two helper functions that we can use to try these ideas out.  \n",
    "\n",
    "The first helper function plots the along-track slopes for each segment:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_segs(D6, ind=None, **kwargs):\n",
    "    \"\"\"\n",
    "    Plot a sloping line for each ATL06 segment\n",
    "    \"\"\"\n",
    "    if ind is None:\n",
    "        ind=np.ones_like(D6['h_li'], dtype=bool)\n",
    "    #define the heights of the segment endpoints.  Leave a row of NaNs so that the endpoints don't get joined\n",
    "    h_ep=np.zeros([3, D6['h_li'][ind].size])+np.NaN\n",
    "    h_ep[0, :]=D6['h_li'][ind]-D6['dh_fit_dx'][ind]*20\n",
    "    h_ep[1, :]=D6['h_li'][ind]+D6['dh_fit_dx'][ind]*20\n",
    "    # define the x coordinates of the segment endpoints\n",
    "    x_ep=np.zeros([3,D6['h_li'][ind].size])+np.NaN\n",
    "    x_ep[0, :]=D6['x_atc'][ind]-20\n",
    "    x_ep[1, :]=D6['x_atc'][ind]+20\n",
    "\n",
    "    plt.plot(x_ep.T.ravel(), h_ep.T.ravel(), **kwargs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The second helper function calculates how close the endpoints of each segment are to centerpoints of that segment's neighbors:\n",
    "\n",
    "<img src=\"images/dh_segment_sm.jpg\"  width=500 height=450>\n",
    "\n",
    "If we look at the difference between $h_m$ - 20 dh_fit_dx and $h_{m-1}$ and at the difference between $h_m$ + 20 dh_fit_dx and $h_{m+1}$, we can get an idea of whether segments $m-1$, $m$, and $m+1$ are consistent.  We'll report the minimum of $|h_m - 20 dh\\_fit\\_dx - h_{m-1}|$ and $|h_m + 20 dh\\_fit\\_dx - h_{m+1}|$.  A small value indicates that a segment is consistent with at least one of its neighbors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "def min_seg_difference(D6):\n",
    "    \"\"\"\n",
    "    seg_difference_filter: Use elevations and slopes to find bad ATL06 segments\n",
    "    \n",
    "    \n",
    "    Inputs: \n",
    "        D6: a granule of ATL06 data, in dictionary format.  Must have entries:\n",
    "            x_atc, h_li, dh_fit_dx\n",
    "        \n",
    "    Returns:\n",
    "        delta_h_seg: the minimum absolute difference between each segment's endpoints and those of its two neighbors\n",
    "    \"\"\"\n",
    "    h_ep=np.zeros([2, D6['h_li'].size])+np.NaN\n",
    "    h_ep[0, :]=D6['h_li']-D6['dh_fit_dx']*20\n",
    "    h_ep[1, :]=D6['h_li']+D6['dh_fit_dx']*20\n",
    "    delta_h_seg=np.zeros_like(D6['h_li'])\n",
    "    delta_h_seg[1:]=np.abs(D6['h_li'][1:]-h_ep[1, :-1])\n",
    "    delta_h_seg[:-1]=np.minimum(delta_h_seg[:-1], np.abs(D6['h_li'][:-1]-h_ep[0, 1:]))\n",
    "    return delta_h_seg\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's see how each of these helps find bad data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "11cc58b46efb40fbb510ff1c5a9264bd",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig=plt.figure(figsize=[6, 6])\n",
    "h1=fig.add_subplot(211)\n",
    "plt.plot(D3['heights']['x_atc'], D3['heights']['h_ph'],'k.', markersize=0.25)\n",
    "# plot the segments:\n",
    "plot_segs(D6, color='b', label='atl06 segments')\n",
    "good=D6['atl06_quality_summary']==0\n",
    "#plot the bad-flagged points\n",
    "plt.plot(D6['x_atc'][~good], D6['h_li'][~good],'y.', label='atl06_quality_summary==1')\n",
    "plt.plot(D6['x_atc'][good],  D6['h_li'][good],'r.', label='atl06_quality_summary==0')\n",
    "plt.legend()\n",
    "plt.title('flagging with atl06_qualiy_summary')\n",
    "\n",
    "h2=fig.add_subplot(212, sharex=h1, sharey=h1)\n",
    "delta_h_seg=min_seg_difference(D6)\n",
    "plt.plot(D3['heights']['x_atc'], D3['heights']['h_ph'],'k.', markersize=0.25, zorder=1)\n",
    "plot_segs(D6, color='b', label='atl06 segments', zorder=2)\n",
    "cm=plt.scatter(D6['x_atc'], D6['h_li'], 4, c=delta_h_seg, vmin=0, vmax=1, cmap='autumn', zorder=3); \n",
    "plt.colorbar(cm, ax=[h1, h2],label='min_seg_difference, m')\n",
    "plt.title('min_seg_difference');\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Selecting segments with atl06_quality_summary==0 or selecting segments with min_seg_difference < 1 is a good way to remove bad points from a plot like this, and combining the two can reduce the blunder rate further.  We'll use the summary flag in this tutorial, but keep both in mind!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. Repeat track data from ATL06\n",
    "Next, let's look at the big picture, combining data from multiple tracks around Pine Island Glacier.  We'll use the projected coordinates and plot on top of the image mosaic.\n",
    "\n",
    "First, we'll every 25th point from one of the center-beam pairs for all files in our data directory.  We'll print an error if the reading fails, but will let the code continue anyway:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "ename": "OSError",
     "evalue": "Unable to open file (unable to open file: name = 'A', errno = 2, error message = 'No such file or directory', flags = 0, o_flags = 0)",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mOSError\u001b[0m                                   Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-40-eee394515130>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m      5\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mfile\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mATL06_files\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      6\u001b[0m     \u001b[1;32mtry\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 7\u001b[1;33m         \u001b[0mD_dict\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mfile\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0matl06_to_dict\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfile\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;34m'/gt2l'\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mindex\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mslice\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m-\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m25\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mepsg\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m3031\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      8\u001b[0m     \u001b[1;32mexcept\u001b[0m \u001b[0mKeyError\u001b[0m \u001b[1;32mas\u001b[0m \u001b[0me\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      9\u001b[0m         \u001b[0mprint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34mf'file {file} encountered error {e}'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m<ipython-input-33-0614904ba59b>\u001b[0m in \u001b[0;36matl06_to_dict\u001b[1;34m(filename, beam, field_dict, index, epsg)\u001b[0m\n\u001b[0;32m     25\u001b[0m     \u001b[0mD\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;33m{\u001b[0m\u001b[1;33m}\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     26\u001b[0m     \u001b[0mfile_re\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mre\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcompile\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'ATL06_(?P<date>\\d+)_(?P<rgt>\\d\\d\\d\\d)(?P<cycle>\\d\\d)(?P<region>\\d\\d)_(?P<release>\\d\\d\\d)_(?P<version>\\d\\d).h5'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 27\u001b[1;33m     \u001b[1;32mwith\u001b[0m \u001b[0mh5py\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mFile\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfilename\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;34m'r'\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;32mas\u001b[0m \u001b[0mh5f\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     28\u001b[0m         \u001b[1;32mfor\u001b[0m \u001b[0mkey\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mfield_dict\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     29\u001b[0m             \u001b[1;32mfor\u001b[0m \u001b[0mds\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mfield_dict\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mkey\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mC:\\Anaconda3\\lib\\site-packages\\h5py\\_hl\\files.py\u001b[0m in \u001b[0;36m__init__\u001b[1;34m(self, name, mode, driver, libver, userblock_size, swmr, rdcc_nslots, rdcc_nbytes, rdcc_w0, track_order, **kwds)\u001b[0m\n\u001b[0;32m    406\u001b[0m                 fid = make_fid(name, mode, userblock_size,\n\u001b[0;32m    407\u001b[0m                                \u001b[0mfapl\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mfcpl\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mmake_fcpl\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mtrack_order\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mtrack_order\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 408\u001b[1;33m                                swmr=swmr)\n\u001b[0m\u001b[0;32m    409\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    410\u001b[0m             \u001b[1;32mif\u001b[0m \u001b[0misinstance\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mlibver\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mtuple\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mC:\\Anaconda3\\lib\\site-packages\\h5py\\_hl\\files.py\u001b[0m in \u001b[0;36mmake_fid\u001b[1;34m(name, mode, userblock_size, fapl, fcpl, swmr)\u001b[0m\n\u001b[0;32m    171\u001b[0m         \u001b[1;32mif\u001b[0m \u001b[0mswmr\u001b[0m \u001b[1;32mand\u001b[0m \u001b[0mswmr_support\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    172\u001b[0m             \u001b[0mflags\u001b[0m \u001b[1;33m|=\u001b[0m \u001b[0mh5f\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mACC_SWMR_READ\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 173\u001b[1;33m         \u001b[0mfid\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mh5f\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mopen\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mname\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mflags\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mfapl\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mfapl\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    174\u001b[0m     \u001b[1;32melif\u001b[0m \u001b[0mmode\u001b[0m \u001b[1;33m==\u001b[0m \u001b[1;34m'r+'\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    175\u001b[0m         \u001b[0mfid\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mh5f\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mopen\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mname\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mh5f\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mACC_RDWR\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mfapl\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mfapl\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mh5py\\_objects.pyx\u001b[0m in \u001b[0;36mh5py._objects.with_phil.wrapper\u001b[1;34m()\u001b[0m\n",
      "\u001b[1;32mh5py\\_objects.pyx\u001b[0m in \u001b[0;36mh5py._objects.with_phil.wrapper\u001b[1;34m()\u001b[0m\n",
      "\u001b[1;32mh5py\\h5f.pyx\u001b[0m in \u001b[0;36mh5py.h5f.open\u001b[1;34m()\u001b[0m\n",
      "\u001b[1;31mOSError\u001b[0m: Unable to open file (unable to open file: name = 'A', errno = 2, error message = 'No such file or directory', flags = 0, o_flags = 0)"
     ]
    }
   ],
   "source": [
    "# find all the files in the directory:\n",
    "ATL06_files='ATL06_20190128163325_04770206_003_01.h5'\n",
    "D_dict={}\n",
    "error_count=0\n",
    "for file in ATL06_files:\n",
    "    try:\n",
    "        D_dict[file]=atl06_to_dict(file, '/gt2l', index=slice(0, -1, 25), epsg=3031)\n",
    "    except KeyError as e:\n",
    "        print(f'file {file} encountered error {e}')\n",
    "        error_count += 1\n",
    "print(f\"read {len(D_dict)} data files of which {error_count} gave errors\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "About 20 of the files had problems, likely because clouds obscured the surface for gt2l.  That's too bad, but we can still work with the data we have. \n",
    "\n",
    "## 3.1. Repeat structure by cycle\n",
    "\n",
    "Let's map the ground tracks for cycles 1 and 2 (not on repeat tracks) and for cycles 3 and later (measured on the repeat tracks)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plt.figure(figsize=[8,8])\n",
    "# hax0=plt.gcf().add_subplot(211, aspect='equal')\n",
    "# MOA.show(ax=hax0, cmap='gray', clim=[14000, 17000]);\n",
    "# hax1=plt.gcf().add_subplot(212, aspect='equal', sharex=hax0, sharey=hax0)\n",
    "# MOA.show(ax=hax1, cmap='gray', clim=[14000, 17000]);\n",
    "# for fname, Di in D_dict.items():\n",
    "#     cycle=Di['cycle']\n",
    "#     if cycle <= 2:\n",
    "#         ax=hax0\n",
    "#     else:\n",
    "#         ax=hax1\n",
    "#     #print(fname)\n",
    "#     #print(f'\\t{rgt}, {cycle}, {region}')\n",
    "#     ax.plot(Di['x'], Di['y'])\n",
    "#     if True:\n",
    "#         try:\n",
    "#             if cycle  < 3:\n",
    "#                 ax.text(Di['x'][0], Di['y'][0], f\"rgt={Di['rgt']}, cyc={cycle}\", clip_on=True)\n",
    "#             elif cycle==3:\n",
    "#                 ax.text(Di['x'][0], Di['y'][0], f\"rgt={Di['rgt']}, cyc={cycle}+\", clip_on=True)\n",
    "#         except IndexError:\n",
    "#             pass\n",
    "# hax0.set_title('cycles 1 and 2');\n",
    "# hax1.set_title('cycle 3+');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that cycles 1 and 2 were not precicely repeated, but all data from cycles 3 and onwards follow an exact set of repeats.  We've labeled the start of each track.  In this area, ascending tracks are labeled on the right (true South) side of the plot, while descending tracks are labeled on the left (true North) side."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.2 A map-view look at elevations\n",
    "Now let's map the elevations associated with these plots."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "# map_fig=plt.figure()\n",
    "# map_ax=map_fig.add_subplot(111)\n",
    "# MOA.show(ax=map_ax, cmap='gray', clim=[14000, 17000])\n",
    "# for fname, Di in D_dict.items():\n",
    "#     # select elevations with good quality_summary\n",
    "#     good=Di['atl06_quality_summary']==0\n",
    "#     ms=map_ax.scatter( Di['x'][good], Di['y'][good],  2, c=Di['h_li'][good], \\\n",
    "#                   vmin=0, vmax=1000, label=fname)\n",
    "# map_ax._aspect='equal'\n",
    "# plt.colorbar(ms, label='elevation');\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Elevations run from low (blue) on the ice shelf, to high (yellow) on the surrounding ridges. There are a few elevation blunders here, but nothing to cry about, and most of the tracks seem to have fairly good coverage."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.3 Plotting a repeat-track elevation profile\n",
    "Next, let's plot a single profile for the black line in the figure above.  We'll plot the surface height (h_li) as a function of the along-track coordinate, x_atc.  Since when we read in the elevations to make the survey plot we only read every 100th elevation, we'll reread the data at full resolution before plotting it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "53e95f6927d44acaa13f0a1aa2efa278",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "No handles with labels found to put in legend.\n"
     ]
    }
   ],
   "source": [
    "D_2l={}\n",
    "D_2r={}\n",
    "\n",
    "\n",
    "# specify the rgt here:\n",
    "rgt=\"0027\"\n",
    "# iterate over the repeat cycles\n",
    "for cycle in ['03','04','05','06','07']:\n",
    "    for filename in glob.glob(os.path.join(data_root, 'PIG_ATL06', f'*ATL06_*_{rgt}{cycle}*_003*.h5')):\n",
    "        try:\n",
    "            # read the left-beam data\n",
    "            D_2l[filename]=atl06_to_dict(filename,'/gt2l', index=None, epsg=3031)\n",
    "            # read the right-beam data\n",
    "            D_2r[filename]=atl06_to_dict(filename,'/gt2r', index=None, epsg=3031)\n",
    "            # plot the locations in the previous plot\n",
    "            map_ax.plot(D_2r[filename]['x'], D_2r[filename]['y'],'k');  \n",
    "            map_ax.plot(D_2l[filename]['x'], D_2l[filename]['y'],'k');\n",
    "        except Exception as e:\n",
    "            print(f'filename={filename}, exception={e}')\n",
    "\n",
    "plt.figure();\n",
    "for filename, Di in D_2l.items():\n",
    "    #Plot only points that have ATL06_quality_summary==0 (good points)\n",
    "    hl=plot_segs(Di, ind=Di['atl06_quality_summary']==0, label=f\"cycle={Di['cycle']}\")\n",
    "    #hl=plt.plot(Di['x_atc'][Di['atl06_quality_summary']==0], Di['h_li'][Di['atl06_quality_summary']==0], '.', label=f\"cycle={Di['cycle']}\")\n",
    "    \n",
    "plt.legend()\n",
    "plt.xlabel('x_atc')\n",
    "plt.ylabel('elevation');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see here that we have elevations from cycles 3, 5, and 6, but cycle 4 is missing.  That's likely because of clouds, and it's a fact of life in places where glaciers are found."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# BREAKOUT: [15 minutes of participant work]\n",
    "\n",
    "We're going to move into breakout rooms for a few minutes.  Introduce yourselves, then take a look at the ATL03 data and at the ATL06 data.\n",
    "\n",
    "1. Take some time to explore the PIG ATL03 data set.  There are two repeat measurements of the same rgt from two different cycles.  We have looked at one beam from each, but you can explore the different beams (strong and weak) and the different pairs (1, 2, and 3.)  Combine the first cell from 2.1.2 with the second cell from 2.1.3 to find examples where:\n",
    "\n",
    "- ATL03 didn't capture the surface\n",
    "- ATL06 is locked on something that's not the surface\n",
    "- The quality-summary filtering strategy didn't work well\n",
    "\n",
    "Please find at least one example of each, and post a screen grab in the #questions channel\n",
    "\n",
    "2. Explore the ATL06 dataset.  The maps in 3.2 show where each track in the data set covers.  You can change the track number in the cell from 3.3, and can use the zoom button to look at different areas of each plot.  Also try looking at different beams (we've mapped only the center pair, but the other two pairs are there as well!)\n",
    "\n",
    "Within your breakout rooms, work together to look at tracks from different parts of the dataset:\n",
    "- inland ice (the sides of the ice stream)\n",
    "- Ice-stream margins with crevasses (found at the edges of the ice stream)\n",
    "- Ridged ice-shelf terrain (roughly the bottom third of the map)\n",
    "- Rifts (cracks within the ice shelf)\n",
    "- Data collected over sea ice\n",
    "\n",
    "Again, post an example to the slack channel, and feel free to post anything else of intrest that you find!\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.4 an repeat-track example with substantial cross-track slope\n",
    "Let's have a look at one of the tracks (rgt 652).  This time we're going to plot both the left and right beams in the central beam pair.  We'll plot the left beam with a dot and the right beam with plus."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "D_2l={}\n",
    "D_2r={}\n",
    "\n",
    "# specify the rgt here:\n",
    "rgt=\"0469\"\n",
    "for cycle in ['03','04','05','06','07']:\n",
    "    for filename in glob.glob(os.path.join(data_root, 'PIG_ATL06', f'*ATL06_*_{rgt}{cycle}*_003*.h5')):\n",
    "        print(filename)\n",
    "        try:\n",
    "            # read the left-beam data\n",
    "            D_2l[filename]=atl06_to_dict(filename,'/gt2l', index=None, epsg=3031)\n",
    "            # read the right-beam data\n",
    "            D_2r[filename]=atl06_to_dict(filename,'/gt2r', index=None, epsg=3031)\n",
    "        except Exception as e:\n",
    "            print(f'filename={filename}, exception={e}')\n",
    "\n",
    "plt.figure(figsize=[5, 8]);\n",
    "ax1=plt.subplot(311)\n",
    "MOA.show(cmap='gray', clim=[14000, 17000])\n",
    "ax2=plt.subplot(312)\n",
    "ax3=plt.subplot(313, sharex=ax2)\n",
    "for filename, Dl in D_2l.items():\n",
    "    print(filename)\n",
    "    good=Dl['atl06_quality_summary']==0\n",
    "    hl=ax2.plot(Dl['x_atc'][good], Dl['h_li'][good],'.', label=f\"cycle={Dl['cycle']}\")\n",
    "    ax1.plot(Dl['x'][::10], Dl['y'][::10], '.', color=hl[0]._color)\n",
    "    Dr=D_2r[filename]\n",
    "    good=Dr['atl06_quality_summary']==0\n",
    "    ax2.plot(Dr['x_atc'][good], Dr['h_li'][good],'+', color=hl[0]._color)\n",
    "    ax1.plot(Dr['x'][::10], Dr['y'][::10], '+', color=hl[0]._color)\n",
    "    \n",
    "    ax3.plot(Dr['x_atc'][::10], Dr['y_atc'][::10],'.', color=hl[0]._color)\n",
    "    ax3.plot(Dl['x_atc'][::10], Dl['y_atc'][::10],'+', color=hl[0]._color)\n",
    "    \n",
    "ax2.legend()\n",
    "for ax in [ax2, ax3]:\n",
    "    ax.set_xlabel('x_atc')    \n",
    "ax2.set_ylabel('elevation');\n",
    "ax3.set_ylabel('y_atc');\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that cycles 3, 4, and 6 are misaligned from one another.  Zooming , we can see that the terrain is higher to the left (plusses) than to the left.  This is very likely biasing the apparent elevation differences.  \n",
    "# 3.5.1 : Correcting repeat elevations for cross-track slope\n",
    "\n",
    "The ATL11 product will soon be available to sort this out, but in the meantime, we can do a basic correction using the across-track slope estimate in ATL06.  Since this is a middle pair, we should be able to use the slope to correct both beams' elevations to y_atc=0.  If it were left or a right pair, we'd correct it back to + or - 3200 m."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig=plt.figure()\n",
    "ax1=fig.add_subplot(211)\n",
    "ax2=fig.add_subplot(212, sharex=ax1, sharey=ax1)\n",
    "color_dict={}\n",
    "for filename, Dl in D_2l.items():\n",
    "    print(filename)\n",
    "    good=Dl['atl06_quality_summary']==0\n",
    "    Dr=D_2r[filename]\n",
    "\n",
    "    hl=ax1.plot(Dl['x_atc'][good], Dl['h_li'][good],'.', markersize=3, label=f\"cycle={Dl['cycle']}\")\n",
    "    # save the color for this cycle so we can use it again\n",
    "    color_dict[int(Dl['cycle'])]=hl[0]._color\n",
    "    atc_corr=Dl['y_atc']*Dl['dh_fit_dy']\n",
    "    ax2.plot(Dl['x_atc'][good], Dl['h_li'][good]-atc_corr[good], '.', markersize=3,color=color_dict[int(Dl['cycle'])])\n",
    "    \n",
    "    good=Dr['atl06_quality_summary']==0\n",
    "    ax1.plot(Dr['x_atc'][good], Dr['h_li'][good], '+', markersize=4, color=color_dict[int(Dl['cycle'])])\n",
    "    atc_corr=Dr['y_atc']*Dr['dh_fit_dy']\n",
    "    ax2.plot(Dr['x_atc'][good], Dr['h_li'][good]-atc_corr[good], '+', markersize=4, color=color_dict[int(Dl['cycle'])])\n",
    " \n",
    "ax1.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The elevations from each cycle now lie directly on top of each other.  Note that there are far fewer elevations from cycle 3, because the dh_fit_dy variable is only available when both beams in a pair have good-quality elevations (based on atl06_quality summary).  \n",
    "\n",
    "### 3.5.2 A look at ATL11 cross-track slope corrections\n",
    "ATL11 takes care of this problem by creating a common surface that can correct the elevations for all available cycles."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "D11=pc.ATL11.data().from_h5(os.path.join(data_root,'PIG_ATL11','ATL11_046910_0306_02_vU07.h5'))\n",
    "fig=plt.figure(); \n",
    "ax_atl11_1=fig.add_subplot(211, sharex=ax1, sharey=ax1)\n",
    "ax_atl11_2=fig.add_subplot(212, sharex=ax1)\n",
    "for col in range(D11.h_corr.shape[1]):\n",
    "    if col+3 in color_dict:\n",
    "        ax_atl11_1.plot(D11.x_atc[:,0], D11.h_corr[:,col],'.', markersize=2, color=color_dict[col+3], label=f\"cycle {col+3}\")\n",
    "        ax_atl11_2.plot(D11.x_atc[:,0], D11.h_corr[:,col]-D11.h_corr[:, 3],'.', markersize=2, color=color_dict[col+3])\n",
    "ax_atl11_1.set_ylabel('height')\n",
    "ax_atl11_1.legend()\n",
    "ax_atl11_2.set_ylabel('height difference from cycle 6')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With a bit of zoom work, you can see that there have been 2 or 3 meters of drawdown between cycle 3 (April-June, 2019) and cycle 6 (October-December, 2020)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
